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Abstract 

We study the bulk properties of isotropic bidisperse granular mix¬ 
tures using discrete element simulations. The focus is on the influence 
of the size (radius) ratio of the two constituents and volume fraction 
on the mixture properties. We show that the effective bulk modulus 
of a dense granular (base) assembly can be enhanced by up to 20% by 
substituting as little as 5% of its volume with smaller sized particles. 
Particles of similar sizes barely affect the macroscopic properties of 
the mixture. On the other extreme, when a huge number of fine par¬ 
ticles are included, most of them lie in the voids of the base material, 
acting as rattlers, leading to an overall weakening effect. In between 
the limits, an optimum size ratio that maximizes the bulk modulus of 
the mixture is found. For loose systems, the bulk modulus decreases 
monotonically with addition of fines regardless of the size ratio. Fi¬ 
nally, we relate the mixture properties to the typical pore size in a 
disordered structure as induced by the combined effect of operating 
volume fraction (consolidation) and size ratio. 

Keywords; Granular mixtures. Discrete element method. Isotropic com¬ 
pression, Extreme size ratio. Rattlers, Elastic moduli 

1 Introduction 

Granular materials are widely used as raw materials in various industries, 
including pharmaceutical, mining, chemical, agricultural, household prod¬ 
ucts and food industries. In many of these applications, processes involving 
milling, segregation, agglomeration, filtration and sieving are common and 
often lead to the generation of granular systems with large size ratios. Deal¬ 
ing with highly polydisperse systems is exceptionally challenging and often 
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requires heuristic assumptions to be made, as prediction/control of the be¬ 
havior is still un unsolved issue. 

On the other hand, it is well known in the geomechanical community 
that the presence of small particles (fines) strongly influences the mechanical 
behavior of granular soils. Scientific work on the topic is extensive. Several 
models have been proposed to describe the variation of stiffness and strength 
of granular mixtures as a function of the volume of fines in the system (see 
[HElEllllISlEllTlEllH] among others). Thevenayagam et al. [I] proposed the 
concept of intergranular void ratio (distinct from the measured, apparent 
voids ratio), assuming that up to a certain fines content (dependent on 
density) the finer grains do not actively participate in the transfer of contact 
frictional forces. 

The problem has been recently approached also numerically [ladiiiisi 
ESI O E], by using the Discrete Element Method (DEM). Ogarko and 
Ending [TT] found numerically that any polydispersity can be replaced by 
an equivalent bidisperse mixture when the size distribution moments are 
matched, in the case of isotropic compression. Recent works by Ueda et al. 
ES] have explored ranges of size and volume ratios of bidisperse granular 
mixtures to evaluate the shear strength in the quasi static regime. Minh 
et al. [13] studied the strong force network in a bimodal granular mixture 
under uniaxial compression, depending on the quantity of fines in the system, 
detecting an optimum value of hne content. Shaebani et al. m used a 
mean held approximation and found a direct relation between the mean 
packing properties of the stiffness components in the case of (uniformly 
distributed) polydisperse aggregates. The micro-macro scaling is realized 
through a combination of moments of the particle size distribution. 

However, all cited works refer to either systems with an homogeneous 
size distribution or to bidisperse mixtures of constant size ratio, where the 
relative volume of particles is varied. To the best of our knowledge, no sys¬ 
tematic study has been done to study the effect of including a small volume 
of hues in a granular aggregate. The interest hrst rises from geophysical 
hazards, like earthquakes, where the material volume remains constant, but 
the size of particles quickly decreases due to breakage. The change in size 
distribution, even limited to a very small volume ratio, have been shown 
to play an important role in soil stability. Applications are uncountable in 
industrial processes, where the focus is optimizing the performances of a 
given (granular) material with minimum modification, i.e., minimum costs. 
Hence, the presence of small particles in a granular mass, which is often 
associated with weakness, is here turned into an asset for functionality. 

We use the DEM to study the effect on micro-and macroscopic quantities 
of a monodisperse granular assembly by substituting only 5/105=4.76% of 
its volume with particles of different size (and same characteristics), thus 
generating a bidisperse mixture, with the main focus on the bulk stiffness. 
We analyze the properties of the granular mixture on two phase spaces: (i) 
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by varying the size ratio of fines to coarse and (ii) by spanning a wide range 
of volume fraction, and find that at each volume fraction corresponds at 
optimum size ratio that maximizes the bulk modulus. 

This paper is organized as follows: The simulation method and parame¬ 
ters used and the averaging definitions for scalar and tensorial quantities are 
given in section The preparation test procedure for creating the granular 
mixture is explained in section]^ Section]^ is devoted to the rattlers (that 
do not contribute to the mechanical stabilities), and the effect of the size 
of fines on microscopic quantities like the coordination number. In section 
we discuss the effect of the size of fines on macroscopic quantities like 
pressure and the jamming volume fraction. The effect of the size of fines 
on the contact network, quantified by the isotropic fabric, is also discussed 
there. Finally, section is devoted to the bulk modulus and its variation 
with the size of fines for different volume fractions. 


2 Numerical simulation and Material properties 


In this section, the procedure of creating the granular mixture is presented. 
Later, we discuss the contact model and the simulation parameters. 

The reference sample consists of = 1050 monodisperse particles A 
with radius = 1.5[mm]. Starting from this base sample, many mixtures 
are created by substituting a given number of — A^J) = 50 particles, 
with particles of species B of different radius vb < va, such that the same 
volume — A^J) (47r/3)r^ = N'^ (47r/3)r^ = of material A is re¬ 
placed by B. The volume ratio of the two components in the final mixture 
is thus: 


4 > = 
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and is much smaller than the pore space of the base material. The size ratio 
is varied systematically from the base case rs/rA = I down to tb/ta = 0.13; 
the number of B particles Nb varies together with r^, while the volume ratio 
is kept constant, as well as the volume of the individual species. The total 
volume of particles is Vr = + = 1.05I/J, so that volume of the box V 

is same for different granular mixtures with different number of B particles 
and at a given volume fraction = VrjV. Note that the substitution 
can be thought of addition when a system containing A^J particles of A is 
mixed with B with volume fraction = 5% of that of A. 

In order to characterize the mixtures with different N'^, we define a 
dimensionless quantity /3 as: 






Nj + Nf 


( 2 ) 
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which is the ratio of small particles B to the total number of particle in 
the system. j3 is the input parameter of the simulation and is systematically 
varied to study its effect on the measured micro-macroscopic quantities. For 
small (3, few big particles B are present in the system, while for large /3, many 
smaller particles B are present. The ratio of N'^ to AiJ in terms of /3 is given 
as 


Nl 1-/3’ 


(3) 


and the size (radius) ratio is 


ta 




1/3 


(4) 


The sample made of only A particles is always used as reference case and 
corresponds to the case vb/va = ^- It provides the minimum /3, /3min = 
$/ (1 -|- <1>) = 0.5/1.05 = 0.0476, and hence the minimum The 

variation of the radius ratio tb/ta is reported in Fig. 1(a) and shows a 
monotonic decrease with /3. 

The Discrete Element Method (DEM) |18j has been used extensively 


to study granular materials in biaxial and triaxial geometries [la [2011211 
|22l|23l|2l|25] under general deformation paths involving advanced contact 
models for fine powders |26l l27j . In this work, however, we restrict ourselves 
to the simplest isotropic deformation test and to the linear contact model 
without any friction between the particles. Since DEM is a standard method, 
only the contact model parameters relevant for our simulation are briefly 
discussed. 

The simplest linear no rmal contact force model when two particles i and j 
interact, as shown in Fig. 1(b), is given as fT = ff^ri = {kijS + '^ij5)h., where 
kij is the contact spring stillness, is the contact viscosity parameter, 
6 is the overlap and <5 is the relative velocity in the normal direction h. 
An artificial background dissipation force, f;, = — proportional to the 
velocity Vj of particle i is added (similarly of particle j), resembling the 
damping due to a background medium, as e.g. a fluid. Note that apart 
than the radius, materials A and B have the same interacting properties, 
i.e. stiffness, viscosity and density, see Table 

For a pair of particles i and j with masses rm and rrij , a typical response 
time is the collision duration fj/ = -n j y/kij jmij — {'^ij l2mij Y, where rriij = 


mimj/{mi + mj) is the reduced mass 


In DEM, the integration time- 


step is chosen to be about 50 times smaller than the shortest time-scale 
tc = min [26]. The parameters used in DEM simulations are presented 
in Table For our system, material B sets the DEM time-step, as vb and 
hence itibb is smallest, leading to the smallest The variation of tc 

with /3 can be seen in Fig. |l(a) tc decreases with increasing /3, meaning 
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Parameter 


Symbol Material A Material B 


Number of Particles N'^ 

Radius r 

Particle density p 

Normal stiffness 
Normal Viscosity 7 

Background viscosity 7 ;, 


Vj = 1000 

ta = 1-5 mm 
PA = 2000 [kg/m^] 
k\ = 5.10® [kg/s^] 
1 [kg/s] 

0.1 [kg/s] 


varied [50-22500] 

VB/rA = 

PB = PA [kg/m®] 

un _ un 

1 [kg/s] 

0.1 [kg/s] 


Table 1: Summary and numerical values of particle parameters used in the 
DEM simulations. /3 is the ratio of particles of B to the total number of 
particles, defined in Eq. Q. $ = 0.05 is the ratio of volume of B to that of 
A in the final mixture. 


that the smaller particles in the mixture lead to a reduction in the collision 
time and hence to a finer time-step. Due to computational limitations, the 
simulations were performed up to /3 = 0.957. 

3 Preparation and test procedure 

Each mixture, made of materials A and B as introduced in section is 
created and further compressed using a unique, well defined protocol. The 
preparation consists of three parts: (i) randomization, (ii) isotropic compres¬ 
sion, and (iii) relaxation, all equally important to achieve the initial mixtures 
for the following analysis. The initial configuration is such that spherical 
particles of particles A and B, are randomly generated in a 3D box without 
gravity, with low volume fraction and rather large random velocities, such 
that they have sufficient space and time to exchange places and to random¬ 
ize themselves. This granular gas is then isotropically compressed, in order 
to approach a direction independent initial configuration with target volume 
fraction vq = 0.64, sightly below the jamming volume fraction, i.e. the tran¬ 
sition point from fluid-like behavior to solid-like behavior [291 Eol EH |32]. 
Isotropic compression is realized by a simultaneous inward movement of all 
the periodic boundaries of the system, with diagonal strain rate tensor 

-10 0 \ 

0 -1 0 , 

0 0 -1 / 

where is the rate amplitude (e^ > 0 in our convention represents compres¬ 
sion) applied to the walls. This is followed by a relaxation period at constant 
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volume fraction to allow the particles to fully dissipate their energy and to 
achieve a static configuration in mechanical equilibrium, indicated by the 
drop in kinetic to potential energy ratio to almost zero. This relaxed state 
is further isotropically compressed until a target maximum volume fraction 
t'max = 0.82 is achieved. The simulations are continued with negative rate 
amplitude in the unloading mode, until the initial t'o = 0.64 is reached. 

For each mixture, configurations at six different v are picked from the 
unloading branch and relaxed, allowing to dissipate the kinetic energy and 
reach unjammed, non-overlapping stable packings. As an example, we 
show in Fig. isotropic samples with (5 = N'^/ (aJ -|- N'^) =0.075, 0.56, 
and 0.96 for loose and dense samples with volume fraction = 0.69 and 0.82 
respectively. 

4 Microscopic Quantities 

In this section, we present the general definitions of averaged microscopic 
parameters including the coordination number and the fraction of rattlers. 

4.1 Mechanically stable system 

In order to properly link the macroscopic load carried by the sample with the 
active microscopic contact network, all particles that do not contribute to 
the force network are excluded from the analysis. Frictionless particles with 
less than 4 contacts are thus ‘rattlers’, since they are not mechanically stable 
and hence do not participate to the force transmission [28l [33l IMj • From 
the snapshots in Fig. where number of contacts of particle p is Cp > 0, all 
the particles with less than 4 contacts are removed. The rattlers exclusion 
is an iterative process until all remaining particles have at least 4 contacts 
{C* > 4), that provides us a completely mechanically stable systems as 
shown in Fig. (only particles B are shown). 

Unless mentioned explicitly, we will denote and Nb as the number 
of mechanically stable particles of A and B, respectively, and use them to 
compute micro- and macroscopic quantities. Any superscript ‘ T’ relates to 
the total number of particles of A and B, including the rattlers. 

4.2 Rattlers 

In Fig. 1^ we plot the number ratio of participating particles B with respect 
to A after removing the rattlers, i.e., Nb/Na- For all cases, the assembly 
contains 95% by volume of big particles of A. Thus, Na after removing 
rattlers is close to AJ, i.e. Na/N'^ ps 1. With decreasing size of B, i.e., 
increasing /3, an initial increase in the ratio Nb/Na is seen, followed by 

^Configurations from the unloading branch are more reliable since it is much less sen¬ 
sitive to the protocol and rate of deformation during preparation | 3,SI I28| . 
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a maximum and a later decrease for all the volume fractions. For small 
/3, few B particles are present with size comparable to A. With increasing 
j3 and for all v, the ratio Nb/Na increases as more B particles of smaller 
size are introduced in the system while the number of A stays constant. 
For a hxed there is an average void size created by A that can be most 
efficiently filled by an optimal (just fitting) E 0 the optimum size ratio tb/ta 
corresponds to a maximum in Nb/Na- Indeed, when /3 increases further, 
meaning more smaller particles B in the system, the number of active (non¬ 
rattlers) particles B decreases, as most of them become rattlers, ‘caged’ in 
the voids of A m- This can be seen comparing Fig. |2(c)| with Fig. |3(c)[ 
Therefore, the ratio Nb/Na decreases after the maximum. 

Another important observation is that with increasing the maximum 
in Nb/Na occurs at increasing (3. This is because for increasing v, the base 
assembly A is more compressed with smaller void size. That is smaller B 
particles already fill efficiently the voids of A as non-rattlers. For the densest 
case u = 0.82, the ratio Nb/Na seems to saturate for large (3. However 
when /3 —)• 1, A^b will decrease and hence the ratio Nb/Na- Due to the 
computational limitations, this observation can not be presented. 

Fig. I^also shows the ratio Nb/Na including rattlers, same for all density 
represented by Eq. Q. N'^/N'^ is higher than Nb/Na, since Nb is smaller 
than N'^, while Na/N'^ ~ 1. Note that the dashed line is closer to the dense 
systems, where the majority of B particles are active particles. 


4.3 Coordination number 

The classical definition of coordination number is C = M/N'^ , where M 
is the total number of contacts and N'^ is the total number of particles 
1331 Eg [II]. In order to quantify the active contact network (excluding 
rattlers), we use the corrected coordination number: C* = M/^/N^^ where 
M 4 is the total number of contacts of the A 4 particles with at least 4 contacts. 


i.e. C* > A (see section 4.1). Note that, after excluding rattlers, the number 
of particles left in the system is {Na + Nb) = N^ < N'^ = (A^J -|- N'^) 

Fig. [^shows the evolution of C* with f3 for six different volume fractions. 
As expected, for a given composition (hxed (3), the total coordination num¬ 
ber of the system increases with volume fraction u as the system becomes 
more dense and particles are both closer and better coordinated. For given 
density u, C* decreases continuously with j3, since the number of non-rattler 
particles N 4 increases faster than the non-rattler contacts M 4 ; i.e., C* de¬ 
creases. At high f3, an increase in C* is seen. This is associated with the 
drop in active particles B, Nb, in other words Nb/Na, as shown in Fig. 

®Big particles A create voids filled by B. Different lattice arrangements of A provide 
the void size such that B touches A particles, giving the size ratio rs/rx for the most 
compact packing. A simple approach to measure this ratio for triangular, tetrahedron, 
square and cubic lattices formed by A particle is shown in appendix [A|. 
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4.4 Dimensionless moments 


The average radius and moments are among the fundamental quantities 
needed to characterize the particle size distribution |14) . Given /(r) as the 
particle radii (size) distribution, f{r)dr is the probability to find the radius 
between r and r + dr, with a normalization condition f{r) = 1 [141133| . 
For a bidisperse distribution /(r) = /^lA (r — ta) + /bA (r — r^), where 
Ja = Na/ {Na + Nb) and fs = Nb/ {Na + Nb) are the number fractions 
of A and B and A(r) is the Dirac-delta function. While = AJ + A^, 
without superscript (T denotes the total number of particles), Na+Nb = A 4 
is the total number of particles in the system with at least four contacts 
{C* > 4, see section [4 t| . The moment is (r"') = f^r^f{r). The 
mean particle radius for a bidisperse distribution is thus (r) = rf{r) = 
Ia^a + Jb^b and the moment is (r”) = Ja^a + Jb^b with /a + /s = 1 - 
Fig. 6 (a) shows the average radius of the system scaled with the radius 
of A; i.e., {r)/rA excluding rattlers. Starting from 1, {r)/rA decreases with 
increasing 13 due to the presence of smaller B particles This decrease is 
faster for higher v and shows an inverse trend with respect to Nb/Na in Fig. 
1^ For /3 —>■ 1, the size of B becomes very small compared to A so that they 
do not contribute, and the system excluding rattlers is mainly composed of 
A, so that {r)/rA -A- 1. 

The dimensionless moments of a polydisperse assembly Oi and O 2 are 
defined as m- 

Oi = VV and 02 = ^ 72 - (5) 

(r- 


’ 


where it was shown that Oi and O 2 are needed to completely quantify the 
fluid-like behavior of a granular assembly well below jamming. For our 
system, Na and Nb change with f3 and volume fraction v, hence Oi and O 2 
are different for different volume fractions. If the dimensionless moments Oi 
and O 2 are known, the 2 “'^ and 3'’’^ dimensionless moments (moment scaled 
by (r)) are: 


^_02 , ^_02 
(r)^ O? (r)3 or 


( 6 ) 


The moment is always greater or equal to the power of the mean 
radius, i.e. (r”)/(r)” > 1 . Therefore, Of < O 2 and Oi < 1, as also shown in 
Ref. [14| . For the full system including the rattlers, 


(r)^ = {l-/3)rA+l3rB and (r^)^ = (l-/3)r^+,dr| and (r^)^ = (l-/3)r^+,dr|. 

(7) 

'‘This is understood from the inequality: {r)/rA = fA + fsrB/rA ~ I — fsA fErs/rA = 

1 — /s(l — rs/rA) < 1, since the second term is positive and smaller than unity. 
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Therefore, Oi and O 2 for the full system become 


Oi = 


[(1 - /5) + /3(rs/rA)] [(1 - ^) + ^{rsIrAf 


213 


(1 - /3) + (5{rBlrAY 


and On = 


Inserting the radius ratio tbI^a = from Eq. (4), where 

0.05 = const, (see section]^, Eq. Q can be re-written as 


[{1- + ^{rB/rA) 

[(l-/3) + /3(rB/rA)3]2' 

( 8 ) 


4> = 


or = 


[(1 _ pfO + $1/3^2/3] [(1 _ ^)l/3 + $2/3^1/3] 
(1 + ^^) 


and O 2 = 


[( 1 _; 3 ) 1/3 + ^ 2 / 3 ^ 1 / 3 ]: 


(9) 

The asymptotic values for and when /3 —)• 1 are <I>/(1 -|- <I>) and 
4>^/(l -|- 4>)^, respectively. Eig. 6(b) and 6(c) show the evolution of Oi and 
O 2 with /3 and u. Both Oi and O 2 are smaller than unity m, decreasing 
with /3, and show a similar trend as {r)/rA- The dashed lines in Eig. 
represent the granular assembly including the rattlers, i.e. only a single line 
for all volume fractions. 


(l + 4>)2 


4.5 Corrected volume fraction 

It is interesting to look closer at the behavior of the corrected volume fraction 
ly*, i.e., the volume fraction excluding the non-active particles 


12 * = N 4 


dvr (r^) 


( 10 ) 


Eig. 7(a) shows the corrected volume fraction v* versus (3. Eor any volume 
fraction v, v* decreases continuously with /3, since the volume fraction of 
rattlers (mainly B) increases with j3. Eor decreasing size of B, more and more 
B particles are ‘caged’ between the big particles A without having sufficient 
{C* > 4) contacts [H]. Eor the reference case, when the radius of B is equal 
to that of A, (leftmost data points), v* k, i/, and for the density close to 
the jamming point, v* ~ 0.98i^, as approximately 2% of the particles are 
rattlers, in agreement with the values reported in ffH for the monodisperse 


case. 


Eor each mixture, we extract the jamming point t'c) be., the volume frac¬ 


tion u when the pressure p of the mixture (defined in section 5.2) approaches 
zero. Eig. |7(b)| shows Vc increasing with /3 and saturating for /3 —)> 1. This 
can be understood since the number of non-rattler particles decreases with 
/3, as also seen in Eig. 7(a), until they reach the number of Therefore, 


with increasing (3, one needs to compress the system further (or increase 
the volume fraction) to make sure particles achieve an overlapping, jammed 
configuration, leading to increase in Uc with [3. Eor /3 —>■ 1 and volume 
fractions near 12 ^, all particles B are rattlers and therefore Vc saturates for 
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/3 —)• 1. Note that v* in Fig. 7(a) [ excludes the rattlers while Vc in Fig. |7(b)| 
includes the rattlers. The relation between Vc and j3 can be fitted by the 
linear relation: 

‘ 7 ^. ( 11 ) 

Pmin 


Vc = T^c+ - ^c) 


where = 0.646 for the monodisperse case for /3min = 0.05/1.05 = 0.0476. 

= 0.682 is the only fit parameter for (3 = 1, obtained by fitting the 
simulation data in Fig. 7(b) up to /? = 0.8. Note that the = 0.646 value 


measured from the simulation for the monodisperse case is consistent with 
results in Ref. m for different system size. Assuming that for /3 —?• 1 only 
A particles contribute to te structure, we estimate the saturation volume 
fraction as 


= 1.0 (1 + $) = 0.678, 


( 12 ) 


in consistent with the measured values for /3 —>■ 1 as shown in Fig. |7(b) 


5 Macroscopic Quantities 

In the previous section, we focused on the averaged microscopic quantities; 
rattlers and coordination number. Next we focus on defining the averaged 
macroscopic quantities - stress and fabric (structure), that reveal interesting 
bulk features and provide information about the state of the packing via its 
response to applied deformation. 


5.1 Fabric 


From the DEM simulations, one can determine the fabric tensor in order to 
characterize the geometry/structure of the static aggregate [35l [36l l33l [T7| . 
defined as 


1 




C„ 


p=i 


C=1 


n 


(13) 


where Vp is the volume of particle p, which lies inside the averaging system 
volume V, and is the normal unit branch vector pointing from center of 
particle p to contact c. Cp is the number of contacts of particle p and N'^ 
represents the total number particles. In the case of isotropically compressed 
systems, the isotropic fabric FJ' is the quantity of interest and is obtained 


by taking the trace of Eq. (13) as: 


Fj = tr (F^) = 


N'^ Cp N'^ 

' ^K.Wl = i 




Tr / . 


(14) 


Note that we exclude iteratively the rattlers from the system (see section 
[^, and observe that their contribution to the fabric is small (as shown in 
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Fig. [8]). Therefore, the isotropic fabric for non-rattler particles with stable 
non-rattler contacts {C* > 4) is: 


N 4 




p=i 


1 

V 


NT 

^pCp, 

p=i 


(15) 


where N 4 = {Na + Nb) < (A^J + N'^) = is the total number of particles 
excluding the rattlers, as defined in section 

Fig. shows the evolution of calculated using Eq. (15) with P for six 
volume fractions. The first important observation if that the contribution 
of rattlers to the fabric is small and is very close to F^. For a given 
mixture (fixed /3), F^ increases with volume fraction meaning that the 
system becomes more connected. On the other hand, for a given v, F^ first 
increases and then decreases, The maximum of TV is correlated with the 
average voids created by particles A, that, for a given i' and size of B can be 
optimally filled (see appendix . The behavior of TV is similar to that of 
number of non-rattlers particles of B, Nb, observed in Fig.|^ since F^ oc N 4 
(Eq. ([^). 

We are interested in the relation of isotropic fabric with the system’s 
mean packing properties e.g. volume fraction, average coordination num¬ 
ber. An expression for isotropic fabric in terms on mean packing properties, 
similar to that given in Refs. [iTiiMKn] , as: 


Fy 


931 ^* 0 *, 


(16) 


where u* and C* are the volume fraction and mean corrected coordination 
number of the system respectively, excluding the rattlers, as defined in sec- 
and <73 is related to the moments of size distribution Eor a bidis- 


tion 


4.1 


perse size distribution, excluding the rattlers is given as (see appendix [B| 
for derivation): 


93 ■■= 93 = 


1 r\n/fA + r%n^^fBX _ {r^)i 






(17) 


where and are the radii of A and B with number fraction Ja and 

Jb respectively. n(r) = 27r 1 — / {f + is the space angle 

covered on a particle of radius r by neighboring particles of radius (r). y 
is the ratio of the linear compacity (or the total fraction of shielded surface 
which is proportional to product of space angle and number of non-rattler 
contacts, defined in appendix [B|) of B to A. The unknown in the functional 


^Refs. |17l IS.ti ITT| used the corrected coordination number C = M 4 /N in Eq. (16 1 , 
while is different than C* = M 4 /W used in this study. It was checked that using C in 


Eq. (16 1 instead of C* worsen the comparison with EV. 
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form of Eq. © is X) the ratio of the linear compacities of B to A (see 
appendix U. Fig.[^ 


shows the evolution of the measured ratio x with 
the size ratio vaItb, as extracted from simulations for different volume 
fractions, x increases with increasing size ratio and is dependent on the 
volume fraction v, in agreement with Ref. [IZIEH]. For fitting the data in 
we propose 


Fig. 9(a) 


1 fl 


1 + erf ( a 




(18) 


where erf(...) is the error function and a{v) = 0.25 {y and b{v) = 4.5 + 
a{v) are empirical relations. 1/A ~ 1/0.4 is the maximum compacity ratio 
max(x) = 1/A and is reached near the jamming transition (see appendix [B| 
for the bounds of linear compacity). 

The measured using Eq. (17) is plotted in Fig. 9(b) g^x is greater 
than 1 for all volume fraction, increasing with /3 to a maximum followed 
by a decrease, similar as Nb/Na shown in Fig. g^ measured assuming 
constant linear compacity [SailTllMl, i.e., gl with X = 1 is also plotted in 
Fig. 9(b) shows similar trend as g^ and is higher. Asymptotic analysis for 
/3 —>■ 1, considering all the particles present in the system and with constant 
linear compacity tells us that gj diverges as (1 — in agreement with 

the behavior shown in Fig. |9(b)[ 

shows the relation of F/ with the mean packing properties via 


Fig. 10(a) 


gs using Eq. (16), and a good agreement is observed with small errors up to 
5% for the highest densities and /3 —>■ 1. Modification of the linear compacity 
helps to improve the relation of FV with the mean packing properties, while 
a constant linear compacity assumption works only for low and up to 
intermediate /3, as seen in Fig. |10(b)[ For dense states and high /3, the 
constant linear compacity assumption can lead to up to 45% error in the 
prediction of F/. This is due to the fact that very small B particles are 
present in large numbers, participating in the contact network, so that the 
assumption of the linear compacity independent with particle radii, is not 
valid anymore. The better understanding of the linear compacity that can 
account for large numbers of very small particles in highly polydisperse 
systems is subject of a future study. Finally, we attribute the poor a greeme nt 
between FJ and gli'C as used in Ref. IMUMUm, as shown in Fig. |l0(c)[ to 


the fact that homogeneous size distributions were used, not excluding one 
of the species strongly from the contact network. 


5.2 Pressure 

Besides the fabric, one can determine the static stress tensor as 

1 ^ 

= (19) 

C=1 
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which is the sum of the dyadic products between the branch vector P = 
and the contact force over all the contacts (an exemplary two 

particle contact is shown in Fig. |l(b)[ ) in the system volume V, where the 
contribution of the kinetic energy has been neglected [35l |28] . The isotropic 
component of the stress is the pressure P = tr((T)/3. The non-dimensional 
pressure is defined as: 


Pn = 


2{rA) 

3k 


tr(a) 


k ’ 


( 20 ) 


scaled by constant 2{rA)/k. In bi-axial experiments, the pressure P can 
be measured, and hence pn can be estimated. Note that pn is used in the 
following to avoid dimensions of pressure. The size sensitive non-dimensional 
pressure is defined as [331 ESI E]: 


2 (r) , , 2 (r) 

P = ^7y^tr((T) = —^P. 


3k 


k 


( 21 ) 


Note that in this work, the pressure calculated considering M 4 non-rattler 
contacts is very close to the one from M contacts, as M — are temporary, 
very weak (rattler) contacts that barely contribute to the average stress. 


Fig. 11(a) shows the evolution of the non-dimensional pressure pn with 
/3 for six different volume fractions u as shown in the legend. For a given 
/3, Pn increases with v as the particles are more compressed [28]. On the 
other hand, for a given density, pn systematically decreases with /3. This 
observation is linked to the behavior of the corrected volume fractions v* 
(section [T^ , as seen in Fig. 7(a), which also decrease systematically with (3. 
For moderate and large /3, most of the contribution to pressure comes from 
particles A, while the contribution of B is negligible, as both overlap and 
radius are small and hence is their stress (proportional to both) becoming 
negligible for large (3 (data not shown). Fig. 11(b) [ shows the evolution of the 
non-dimensional pressure p using Eq. (m m- For the smallest {3, where 
the system is composed of only A particles pn and p are the same. For fixed 
v, p decreases much faster than pn with /3, since p is a product of pn and 
(r), both decreasing with (3. In this case, the behavior of p differs from TV, 
see Fig. [^ Another important observation is that the behavior of p has a 
similar trend as C* in Fig. For a given mixture (fixed /3), TV increases 
with volume fraction i/, as do both coordination number and pressure, 

We try to better understand the evolution of the non-dimensional pres¬ 
sure by looking at the individual components that contribute to Eq. ( [19| ). 
Due to the linear contact model used without any tangential component, 
the force and the branch vectors are parallel for all contacts (see Fig. [1(b)]). 
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Hence, Eq. (21) becomes 133 E]: 


p = 


2 (r) 

3k 


tr(o-) 



2 (r) 

3H 


Mi 

^ tr (/cn'" 

C=1 




2 / \ -^4 

= ^ E [('=) + 'y [<*-> + <1 ?■ [('=> 

C=1 


3E 


2 ^ 

3E 
C*v* {r\ 


M4 [(U(5c) + 

,3^ + {I'A)] 


=1 


_(47r/3)iV4(r3)/E_ 


( 22 ) 


27r (r' ^ 

where the rattlers offer no contribution and the prime ' represents the fluc¬ 
tuations with respect to the average. The first term in Eq. (22) considers 
the average overlap ((5c) and the average branch vector (/c). The second 
term is the contribution due to the correlated fluctuations in branch vector 
and overlap. 


Fig. 11(c) shows the evolution of the first term {r){lc){Sc)/ A) in Eq. 
). For a given /? with increasing u, the average overlap (Sc) increases |11| . 
(Ic) slightly decreases and (r)/(r^) increases. Therefore {r){lc){Sc)/ A) in¬ 
creases with L' (for fixed /3) and decreases with (3 (for fixed v), as seen in Fig. 
except for very high z/’s and /3’s. The fluctuation factor {r){l'A)/ A) 

^ The common term C* 

has a similar trend as of (r) {Ic) {6c) / {r^'' 


11 c 


increases with both u and /3, as seen in Fig. 11(d) 

seen in Fig. Comparing Fig. 


as 


11 c 


and Fig. |ll(d)[ we conclude that the decrease in p with (3 observed 
in Fig. 11(a) is mainly associated with the decrease of both {r){lc){6c)/(A 
and C*, while the fluctuation term is very small. 

The non-dimensional pressure p can be written in the same form as given 
in Ref. [33] : 

p = Po^^-^(-ev) [1-7p(-ev)] (23) 

where the quantity (—ffv) is the true or logarithmic volume change of the sys¬ 
tem proportional to the ratio of average overlap to the mean radius {6c) / (r) 
(see appendix [C|) . po ^ 0.043 for uniform size distributions 1371 ITTl 135] : 
however, po is not constant for an arbitrary wide bidisperse size distribu¬ 
tions, as the case in this work, as shown in Fig. 12(a) (see appendix [C| for 
more details about calculating pq). Fig. 11(b) shows the prediction for the 
non-dimensional pressure p using Eq. (|23[) without the second small term. 
Since, 7 p is positive, the pressure is slightly over predicted, mainly in the 
dense regime and for the monodisperse case. Finally, Fig. |12(b)| shows a 
perfect prediction of the measured pressure p when compared with Eq. (23), 
again without the non-linear term. Only for highly dense cases, a maximum 
error of few percent can be seen, which could be avoided by including the 
non-linearity with 7 p. 
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As final stage, we want to study how the stiffness of the granular assembly 
varies with the contribution of the fines. 


6 Bulk modulus 

For each granular mixture, examples are displayed in Fig. we calculate 
the bulk modulus by first relaxing (see section and then applying an 
incremental pure volumetric perturbation of small amplitude to the sam¬ 
ple {diy ^ 0.00015) [37j. The bulk modulus is then the ratio between the 
measured change in pressure and the applied strain dv/v, small enough to 
prevent irreversible contact rearrangements m- 


K' = 1^ 


cLP 


The non-dimensional bulk modulus is thus ESEZIESI: 

K = 

k du 


(24) 


(25) 


where 2(rA) is the average particle radius of A which provides the backbone 
to the granular assembly and k is the particle stiffness. Just like the pressure 
P, the bulk modulus K' can be estimated in the bi-axial experiments, and 
hence K can be calculated. 


Fig. 13(a) shows the evolution of the bulk modulus K plotted against 


/3 for different volume fractions u. As expected, K increases systematically 
with density z^. For loose states {u = 0.69, 0.72), K mostly decreases with 
increasing /3 to the limit case /3 —)> 1, as discussed below. The behavior 
associated with denser states is much more interesting, as we observe an 
increase in AT to a maximum, followed by a decrease for larger /3. Note 
that the value of f3 where K becomes maximum increases with increasing 
densities and the maximum also becomes stronger. From Fig. 13(a) 


we 


extract very important insights: (i) The stiffness of a granular assembly can 
be manipulated by only substituting with a small amount of fines to the base 
material (in this case 5/105%). (ii) We can control the “direction” of the 
change (stiffness enhanced or lowered) and the magnitude of change through 
the density and the size of the small particles, (iii) For loose material, there 
is no enhancement, (iv) For dense material, for a given density, there is an 
ideal size of fines that leads to the maximum in the bulk stiffness. 

We associate the different trends observed for loose and dense systems 
with the ability of the fines (material B) to fill the voids formed by particles 
A. In the loose state, particles B are smaller than the average void size, so 
they act as rattlers, and do not contribute to the force network, leading to a 
decreasing bulk modulus with /3. With increasing density, the void size gets 
smaller and compatible with the size of particles B, and thus contribute to 
the active contact network (see appendix [A]) . 
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The maximum in K observed in Fig. 13(a) is different from that observed 
in the isotropic fabric shown in Fig. Thus in the case of a bidisperse 
granular mixture with wide size ratio, not only the contact network controls 
the bulk modulus, as the case for uniformly polydisperse systems, where the 
bulk modulus is directly associated with HU. 

Fig. 13(b) shows the value of K against v for two extremes: smallest 
/d = /3min and /3 = 1 as shown in the legend, both represent monodisperse 
cases of only A particles. K{I3 = /3min) is the left most data points in 
Fig. 13(a), increasing with u. For /3 = 1, particles B are inhnitely small 
and therefore do not participate in the contact network. Thus, the value 
A(/3 = 1) is obtained by removing B from the system. This leads to a 
slightly smaller volume fraction (1/1.05 times the original) and K{I3 = 1) is 
smaller than K{/3 = /Smin) and the two lines being parallel. The two data¬ 
sets for the limit cases in Fig. 13(b)|show good agreement when fitted by Eq. 


(25). As already visible in Fig. 13(a), Fig 


stiffness measured for each volume fraction from Fig. 13(a 


13(b) shows that the maximum 
does not lie in 


between of the two extremes K{f5 = 1) is smaller than K{f5 = /3min)- 


In order to understand the behavior of K observed in Fig. 13(a) 


we 


look at the contributions to the bulk modulus of the three types of contacts 
present in the system, namely AA, AB and BB. It is straightforward to show 
that K = Kaa + Kab + Kbb- Since the change in stress on B particles is 
very small in average, A'ss is negligible (data not shown) when the granular 
assembly is subjected to small perturbation dv, and hence K ss Kaa + Kab- 
Fig. 14(a)| shows the bulk modulus for AA and AB interactions, Kaa and 


Kab- Kaa remains almost constant with /3, except for the smallest i^, where 
it slightly decreases. Kab remains small for loose system, as B particles 
mostly are rattlers. For high density, we observe an increase in Kab with /I 
followed by a decrease. This signifies that the trend observed for K in Fig. 


13(a) is mainly related to the behavior of AB interactions, while the actual 


value depends on the contributions of AA main network. 

The radius of B, governed by /3 at a particular volume fraction v, plays an 
important role not only in hlling the voids of A, but also in contributing to 
the strong force network, leading to the maxima in bulk modulus K. Now we 
want to relate the bulk stiffness K with the packing properties, in a similar 
fashion of Eqs. (16) and (23) as adopted for fabric and pressure respectively. 
We use the relation proposed in HU to link K to the polydispersity and the 
mean packing properties of the sample [33 [33]: 


2{rA)p*o9si^*C* 

Vc\fW) 


1 27 p ( Cv) -|- ( Cv) (1 7p ( ^v)) 


(91nFv 

a(-ev)_ 


(26) 


where = 0.043, 'jp = 0.2 are constant fit parameters taken from Refs. |28l 
Esmu and the Uc is the jamming volume fraction, gs is the size distribution 
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factor and for our bidisperse distribution is given by (see appendix |T7] : 
(r) + _ {r){r^)g 


Qs = 


{r^) 




{r^) 


(27) 


where the same modification for y as given in Eq. (18) has been adopted 
Fig. 14(b) shows the variation in gg with /3 for different volume fractions 
Qs starts from 1 and decreases with /3, with gg = \ recovered only for 
small densities at larger /3, when the system behaves like an assembly of 
only A particles and B do not contribute to the contact network (rattlers), 
see also in Fig. For dense states, gg decreases continuously and reaches 
0.75. The asymptotic value of gg where also rattlers are considered diverges 
with 1 — /? with power law -1/3, dotted lines in Fig. 14(b)l It is worthwhile 


to notice that the gg in Eq. (27) is different from g^ as used in Eq. ( |17[ ) for 
isotropic fabric Fy, meaning that bulk modulus and fabric depend on the 
polydispersity in a different fashion. In Fig. 13(a)[ the prediction of Eq. 
(26) is reported and an excellent agreement is found. It is noteworthy that 


neglecting the particle polydispersity via gg over-predicts the bulk modulus 


as much as 25%. Note that in Eq. (26) a constant Pq is used, while this is not 
the case for extreme polydispersity in our system (see Fig. 12(a)). Future 
research will focus developing analytical relation for K from Eq. (23), where 
the dependence of po on volume fraction u is also considered. 


7 Qualitative Results 

This paragraph is devoted to the qualitative discussion of the findings of 
this study. 

Starting from a base assembly consisting of particles A, a certain vol¬ 
ume is substituted with smaller size particles B, while the total volume 
is kept constant. We restrict ourselves to a small amount of additives 
(5/105=4.76%), i.e., much less material than would be necessary to fill the 
pore-space in the base material and focus on the effect of particle size-ratio. 
We study the two limits of either similar sizes {vb ~ ta) and of very small 
sized B (r^ <C va), as well as the interesting regime in between the limits. 

Substituting A with similar sized particles B is unlikely to change the 
system properties significantly, since the new particles fully participate in 
the mixture (besides a few rattlers). Thus, we observe a small effect of the 
polydispersity due to the new species on the bulk properties. On the con¬ 
trary, when substituting with very small particles B, the mechanical prop¬ 
erties of the mixture are practically the same as of the base material A 
alone, since the B particles are so small that they can move in between par¬ 
ticles A, freely passing through the pore-throats and thus escaping the pores 

®For a uniform distribution, is very close to unity and is independent of the width 
of the distribution. That may be the reason it did not appear in Refs. |371133| . 
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whenever necessary to reduce their stress. In the intermediate size-regime 
(roughly 1/2 > rB/fA > 1/5) a little volume of particles B can change the 
mixture properties considerably, providing systems with higher mechanical 
bulk modulus as compared to the mere interpolation between the two limit 
cases. 

Assume a pore-size distribution with most pores (formed by A) between 
a cubic and a hexagonal local structure, such that they can accommodate 
particles of sizes between rs = (\/3 — l)rA ~ 0.732r^ and rg = (y^ 3/2 — 
l)r^ Ri 0.225r^, respectively (see appendix for discussion on the pore 
sizes corresponding to different packing arrangements). Thus, the typical 
pore-size is r^/r^ (rg -|- r^j^VA ~ 0.48, while there are no pores outside 
of the range. 

There are two possibilities to scan this range of pores for a given volume 
of substitution. One can either change (i) the size of particle B, rs (or /3) 
or (ii) the size of the pores (by changing the volume of the sample). 

(i) Assuming hxed volume of the system, for rs > rg, any particle B 
sitting in a pore between particles A will mechanically contribute to the 
packing, and we are in the large size of particle B limit. When decreasing 
rs below rg, more and more pores will lose the mechanical contact with 
their (caged or trapped) particles B, until, for vb ~ rg, practically all pores 
filled with single particles B have lost contact with their cages. However, 
pores hlled with more than one particle B still could contribute to the force 
network, so that the number of particles B becomes important (which is 
considerably increasing with decreasing size vb)- When < r 4 = (-v/2 — 
l)rA ~ 0.414ryi also multiple B particles lose their efficiency, since they can 
escape through square pores and for r^ < rg = (2/\/3 — l)rA ~ 0.155r^ 
even through the smallest triangle pores, i.e. we are in the small /3 limit. 

{a) For a fixed size of particles B, increasing the volume fraction (de¬ 
creasing the volume) will reduce the available pore sizes and thus shift the 
whole phenomenology towards smaller tb- All pores become smaller and, for 
the largest densities used, the smallest pore-throats rg are almost closed so 
that the escape mechanism is hindered, but not blocked since there are still 
other (e.g. square) shaped throats in the (disordered, non-crystalline) pack¬ 
ings. The reduction in pore size is proportional to the typical AA-contact 
deformation, which in turn is proportional to the pressure (first order) and 
thus to the density relative to the jamming density. 

Combining the two cases, the maximal bulk stiffness K of the packings, 
as function of volume fraction, occurs at oc rp(l — SaaI^a) ~ ?’p(l — 
hi(i//i/c)). Furthermore, we attribute the increase of the maximal K with 
increasing compression to the reduced mobility of the small particles due to 
smaller pores and pore throats. 

Fig. 15 shows the variation in the bulk modulus K, against tbIta for 
different volume fractions. For dense assemblies, K attains a maxima near 
the size ratio corresponding to tetrahedron lattice, meaning that the dense 


18 



state is more likely to create such a configuration and this is the efficient 
arrangement. For the system becoming looser, the maximum of K moves 
towards high size ratios tbIva corresponding to a cubic-like configuration. 

8 Conclusion and Outlook 

In this study, we use DEM simulations to study the bulk properties of a gran¬ 
ular assembly, initially composed of monodisperse particles. A fixed volume 
of 5% of these monodisperse particles are substituted with 5/105=4.76% of 
fines, in order to create a bidisperse mixture. The focus is on the manipula¬ 
tion of the properties of industrial mixtures with minimal costs/alterations 
by introducing a minimal amount of fines in the assembly. The system is 
characterized by the number ratio /? = fines to the to¬ 

tal number of particles and we study the combined effects of /3 and volume 
fraction on the micro- and macroscopic properties of the mixture. 

Important highlights are extracted regarding microscopic and macro¬ 
scopic (bulk) information of granular mixtures. The static pressure due to 
particle interactions and the coordination number (excluding rattlers) de¬ 
crease monotonically with /I, with small variations for /3 —)> 1. The isotropic 
fabric 7%, a measure of the contact network density, decreases with (5 for 
loose systems (in agreement with pressure behavior), since large pores cre¬ 
ated by big particles provide space for fines to be ‘caged’. The behavior for 
higher densities is different, as TV first increases with /3 and then decreases 
for /I —)• 1. In the first stage, the system is more coordinated and fines effi¬ 
ciently pack the voids, while when /3 —)> 1, most hnes are rattlers, and thus 
Fv decreases. The fabric is well described by the relation Fy = g 3 h'*C*, as 
introduced in Refs. ISgiMlttZlCI!, when 53 is properly modified with a non 
constant compacity accounting for large polydispersity. 

Finally, we focus on the effective bulk modulus K, measured by applying 
small volumetric perturbations to the system. The behavior of K is different 
from both pressure and fabric. For loose systems, a monotonous decrease is 
observed, while for denser systems, K first increases, reaching a maximum, 
whose value depend on the density of the sample, and later decreases. /3 = 1 
can be thought as the case of infinitely small fines and thus resembles the 
monodisperse case with volume fraction 1/1.05 with respect to the original 
case. 

In this study, we focus on the bulk properties of a granular assembly 
when only a small volume fines is included. Next natural step is studying 
the influence of the volume of fines on the material behavior, that is to 
explore different regimes of coarse-fine mixtures. Finally, the focus can be 
moved towards other loading paths (uniaxial compression or shear) to study 
the effect fine content on the evolution of the full elastic tensor. 
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A Radius ratio in different lattice configurations 


In this appendix, we focus on the possible arrangements of particles A and 
B in the granular assembly in order to characterize some special sizes of 
voids in the sample. At a given density, particles A create voids whose 
size depends on the geometry. Among the many possible arrangements of 
A, few possibilities are triangular, tetrahedron, square and cubic lattices as 
shown in Fig. A.l Particles B can efficiently fill these voids at special ratios 


between the void size and the radius of A. 


Fig. A.l.(a) shows a configuration where three A particles are arranged 


in a plane forming a triangular lattice and one particle B of radius is 
located in the void, centered at O, and touches A particles. Thus, 


that means: 


AP = AOcos(ZPAO) = AOcos(30°), 


ta = {rA + rz) 


From here we can obtain the size ratio such that B efficiently fills the pore 
throat formed by three A: 

^ ^ (A.l) 


n 

TA 


— = ^ - 1 « 0.155. 


%/3 


Fig. A.l.(b) shows a sample configuration where A particles are arranged 


in a tetrahedron lattice i.e., a local hexagonal structure, with three of them 
on a plane while the fourth is out of the plane. Assuming the tetrahedron 
is centered at the origin O, that is also the center of particle B of radius 
rg, while AB is one side of the tetrahedron and connects the centers of 
two A particles, {AB = 2rA)- The four vertices of the tetrahedron are 
(±1,0, —1/\/^) and (O, ±1, l/-v/2). The tetrahedral angle ZAOB is 
arccos (—1/3) ss 109.47°. Using the law of cosines we get 

AB^ = AO^ + - 2.AO.ORcos(ZAOB), 

that in terms of radii becomes: 

(rA + TAf = {ta + + {rA + - 2 (r^ + rg) (r^ + rg) cos(arccos (-1/3)) 


= {ta + r^Y 
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The tetrahedron void ratio is thus: 


^ - 1 « 0.225. 


(A.2) 


n. 

rA 


Next in Fig. A.l.(c), we show a configuration where A particles are 
arranged in a planar square lattice and particle B of radius sits in the 
void of A. Using Pythagoras’ theorem, the relation between the sides of the 
lattice is: 

AB'^ = AO^ + OB"^ = 2AO^, 


and introducing the radii 


{rA + rAf = 2{rA + r^f 


and the size ratio for efficient packing is 

^ = V 2 -I. 

ta 


0.414. 


(A.3) 


Finally, Fig. A.l.(d) shows a configuration where A particles are arranged 
in a body centered lattice with particle B of radius rg in the center of the 
cube touching A. Using again Pythagoras theorem, we can write 

= AD^ + CD^ + CG^ = 3AD^, 


and 

{2rA + 2r8)^ = 3 (r^ + , 

that leads to the cubic void size ratio: 

^ = V3 - 1 PS 0.732. (A.4) 

ta 

By comparing the four cases considered here, the triangular lattice pro¬ 
duces the smallest size ratio r'^jrA ~ 0.155 while the cubic lattice r^jrA ~ 
0.732 creates the largest one. 

For the case of overlapping spheres, the size ratio must be corrected by 
including the average overlap between AA ((h^^)) and AB {{Sab)) interac¬ 
tions. 


B Measurement of and gg for fabric and bulk 
modulus 


We are interested in relating the isotropic fabric with the mean packing 
properties (coordination number and volume fraction) via TV = gsi'G. The 


continuous limit of Eq. (15) is given by m- 


JV 

Tv = -^ / V{r)G{r)f{r)dr = g^v*C* 
r JQ 
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(B.l) 





where C{r) is the coordination number of a particle with radius r, volume 
V{r) = 47 rr ^/3 and /(r) is the particle size (radii) distribution defined in 
section [4^ The corrected volume fraction in the continuous form is given 
as; 


N4frVir)f{r) 

V 


(B.2) 


Let’s assume that a reference p-particle with radius r in the system has a 
contact with a neighboring particle of average radius (r). The space angle 
covered by the neighboring particle on the reference particle in a three- 
dimensional packing of sphere is given as [IZII33]: 


n(r) = 27r 



(B.3) 


and the linear compacity (or the fraction of the shielded surface) associated 
with a single interaction is: 

Cs(r-) = (B.4) 

The total compacity of the reference p-particle interacting with its C{r) 
non-rattler neighboring particle of average radius (r) thus becomes: 


Cs{r) = 


dvrr^ 


C(r) 

^Q(r) 

p=i 


r = 


n{r)C{r) 


dvr 


(B.5) 


Cs{r) decreases with r, starting from 1 when r/{r) —>■ 0 and reaches a con¬ 
stant value for r/(r) > 1 . Refs. [T71155] have shown that Cs(r’) decreases 
with increasing particle radii and saturates to a constant value G [ 0 , 1 ] for 
large sized particles. It is also dependent on the volume fraction of the sys¬ 
tem. Large differences in particle number and size ratio affects the linear 
compacity Cs(r). Cs(r) has two bounds: 

i) Upper bound: the maximum compacity is reached when a small par¬ 
ticle is surrounded by two big particles. Therefore, for the lower bound, 
r/(r) —)• 0 leading to Vt{r) = 27r and hence max [cs(r)] = 1 . 

ii) Lower bound: Near the jamming transition (loose states), mainly the 

big particles remain in the system while the smaller particles act as rattlers. 
To be mechanically stable, big particles need six contacts. Using r/(r) —>• 1, 
we have II(r) = 27r (l — \/3/2) and hence min [cs(r)] = A = = 

= 3 (l — \/3/2) 0.4, generally reached by big particles at low 

volume fraction, i.e., near the jamming transition. Using the definition for 
the average coordination number C* in the continuous limit: 


C* = 


r 

C{r)f{r)dr = dvr / 

Jo 


Csir)f{r) 

Q{r) 


(B. 6 ) 
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and using Eq. (B.5), we get 


C* =4tt 


Csir)fir) 

n{r) ' 


Combining Eqs. (B.2), (B.5) and (B.7) in Eq. (B.l), we have: 

r^Cs{r)n{ry^f{r)dr 


93 = 


(B.7) 


(B.8) 


Cs{r)n{r) ^f{r)drf^ r^f{r)dr' 

In a similar fashion, we use a correction term gg as proposed in Ref. m 
to link the bulk modulus R of a granular mixture to the polydispersity as: 

rf{r)dr r‘^Cs{r)Q{r)~^ f{r)dr 


/o°° Cs(r)fl(r) ^f{r)drf^r^f{r)dr 

C Analytical expression for pressure 

In order to better understand the final analytical expressions, the stress is 
rewritten and re-phrased, starting from the traditional definitions. Revisit¬ 
ing Eq. (21), we have: 


(B.9) 
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p=i 


C=1 
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where subscript p and c stand for particles and contacts respectively. 
Average overlap per contact is: 


E Ni ^Cp r V^JV4 e 

_ P=1 2-JC=l _ Z^p=l Z^C=1 "c 

{9c) — 


yNi ^Cp 


Mi NiC* 

Similarly, for the average squared overlap, one can write: 

E Ni r2 \^Ni ^Cp r2 

_ p=l Mc=l _ Mp=l Mc=l 

~ Mi ~ NiC* 

Introducing the average overlap for particle p as: 
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(C.4) 


23 






















where { 6 c) is the average overlap per contact. Eq. (C.l) then can be written 
as: 


P = 
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2 (r) 
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A^4 Cp 
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Ni 
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3E 
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3E 
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(C.5) 


Introducing the normalized particle radius ^p = r^/(r) and overlap Ac = 
6 c/{r) leads to: 

C*u* (Ac) 


P = 


AtT (^3^ 


2 {^p(pp) 


(A|)\ 


(Ac) 7 


C*u* 


where 


9p — 


47r 

{^p4’p 


(Ac) { 2 gp bp{A.c)) 

1 (A^) 


ibp — 


(C. 6 ) 


(C.7) 


(e|) ’" (Cp^) (Ac)2 

The normalized average overlap (Ac) is logarithmically related to the volume 
fraction of the present state via as also presented in Refs. 


(Ac) = Di-e^) = Dln{- 


(C. 8 ) 


Fig. C.l.(a) shows the measured average overlap per contact (A)c against 
ln(i//t'c) and the slope of the linear line is D = 0.425, in consistent with 
the measured value in Ref. [2H]. Therefore, Eq. (C. 6 ) can be written in the 
same form as given in Ref. [33]: 

u*C* 


P = Po -(-ev) [1 - 7p(-ev)] 


(C.9) 


where po = VcgpD/2TT and 7 p = bD/2gp. The unknowns are gp and bp. 
Assuming that total force on a particle is proportional to square of contacts it 
has with neighbors oc Cp, hence c/p Cp/{Cp), where C{v) dvr 

(see appendix]^. Therefore, for a continuous distribution, gp is given as: 

( 0 ^ /o^ rcs{rf9.{r)~‘^f{r)dr 


9p — 


(^^) Cs(r)^fi(r) '^f{r)dr' 


(C.IO) 
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which for a bidisperse size distribution is: 


9p — 


{ry rA^A fA + fe 

(r3) + 


(C.ll) 


where and vb are the radius of A and B with number fraction Ja and /b 


respectively. Note that the second term in Eq. (C.6) is very small (maximum 
10 % for dense volume fractions) and is a subject of future research. Qp 
measured from the simulations is plotted in Fig. C.l.(b)| against Eq. (C.ll) 
and the results are in close agreement. 
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Figure 1: (a) Variation of the radius ratio and smallest collision 

duration tc with /3 = N '^/ (Vj + N'^) . (b) Sketch of two particles in contact 
and the direction of the force and branch vectors. 
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(a) 


(b) 



(e) (f) 

Figure 2: Snapshots of the composite material. White and black particles 
are particles of A (large) and B (small), respectively. Different rows repre¬ 
sent different (3 = N'^/ (A^J -|- = 0.075, 0.56, and 0.96, with size ratio 

fB/rA = 0.5, 0.27 and 0.14. Left and right columns correspond to total 
volume fractions v = 0.69 (loose) and 0.82 (dense), respectively. 
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(e) (f) 

Figure 3: Snapshots of the fines (same as Fig. without rattlers, i.e., 
for clarity only particles B are shown. Different rows represent different 
/3 = N'^/ {N\ + iV^) = 0.075, 0.56, and 0.96, with size ratio vb/ta = 0.5, 
0.27 and 0.14. Left and right columns correspond to total volume fractions 
v = 0.69 and 0.82, respectively. 
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Figure 4: Ratio of B particles with respect A: Nb/Na, when the assembly 
contains no rattlers; the dashed line considers also rattlers (see Eq. (|)),i.e., 
N'^/N'^. All are plotted against the number fraction f3 = N'^/ (iVj + N'^). 
Different colors represent the volume fraction v as shown in the legend and 
the arrow indicates increasing u. 



P 


Figure 5: Coordination number excluding rattlers plotted against the num¬ 
ber fraction /3 = N '^/ • Different colors represent the volume 

fraction u as shown in the legend and the arrow indicates increasing i'. 





(a) 


(b) 


(c) 


Figure 6: (a) Average radius (r) scaled dimensionless moments (b) Oi 
and (c) O 2 , measured using Eq. (§, excluding rattlers, plotted against the 
number fraction /? = N'^ / (A^J -|- A^^) • Different colors represent the volume 
fraction v as shown in the legend. The dashed lines that consider also rattlers 
are Eqs. Q and ([^. The arrows indicate increasing u. 
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p p 

(a) (b) 


Figure 7: (a) Volume fraction v* of the system excluding rattlers plotted 
against the number fraction /? = N'^ / (Vj + . Different colors repre¬ 

sent the volume fraction u as shown in the legend and the arrow indicates 
increasing (b) Evolution of jamming point t'c with (5. The solid line is a 
linear fit to the simulation data using Eq. ©• The dashed horizontal lines 
are: for the smallest /3, fit parameter obtained for /3 —?■ 1 and is 

the estimated value for /3 —)• 1 using Eq. ( |12[ ) when the system contains only 
particles A. 
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Figure 8: Isotropic fabric plotted against the number fraction /3 = 
N '^/ (-^A + ^b') ■ Flifferent colors represent the volume fraction v as shown 
in the legend and arrow indicates increasing v. Open symbols are corre¬ 
sponding FJ that includes the rattlers as well. 
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(b) 


Figure 9: (a) Ratio of linear compacities of particles B to A, x, versus the 
radius ratio r^lrB measured from the DEM simulations (symbols) and the 
corresponding line is the analytical fit to the data using Eq. (18). The dashed 
line is constant linear compacity, i.e., x = 1- (b) for fabric calculated 
using Eqs. 0 with X estimated using (solid symbols) and measured 
from constant linear compacity assumption x = 1 (small open symbols), 
plotted against the number fraction /3 = N '^/ (AJ + A^) excluding the 
rattlers. Different colors represent the volume fraction i' as shown in the 
legend. The dashed line is g^ and considers also rattlers. 
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Figure 10: Isotropic fabric plotted against Eq. (16) with calculated 
using Eq. 0 with (a) non-constant linear compacity, with y computed 
based on volume fraction and radius ratio using Eq. (18) and (b) constant 
linear compacity, with y = 1. The dashed black line has slope 1. Different 
colors represent the volume fraction i' as shown in the legend, (c) The total 
isotropic fabric including the rattlers EJ compared to the relation presented 
in Ref. |33j . and the arrow indicates increasing /3. 
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0.691 0.719 0.747 0.772 0.794 0.820 T- 0.691 0.719 0.747 0.772 0.794 0.820 T- 



Figure 11: (a) Non-dimensional pressure (scaled by constant 2{rA)/k) pn, (b) 
non-dimensional pressure p scaled by 2{r)/k (solid symbols) and prediction 
using Eq. (23) (open symbols) (c) product of mean radius (r), branch vector 
{Ic) and particle overlap {6c) scaled with the third moment (r^) and (d) 
product of mean radius (r) and the corresponding fluctuation term 
scaled with the third moment (r^), calculated using Eq. (22), plotted against 
the number fraction /3 = N '^/ (-/Vj -|- N'^'j. Different colors represent the 
volume fraction i' as shown in the legend and arrows indicate increasing ly. 
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Figure 12: (a) Measured po plotted against the number fraction (3 = 

N '^/ {^A + ^b) • dashed line represent the constant value of pq for 

the monodisperse case, (b) Non-dimensional pressure p plotted against Eq. 
( |23[ ) without the second term. The dashed line is the linear line with slope 1 
and the arrow indicates increasing ly. Different colors represent the volume 
fraction ly as shown in the legend. The arrow in (b) represents increasing (3. 



(a) (b) 


Figure 13: (a) Bulk modulus (scaled by constant 2{rA)/k) K measured using 
Eq. (25) (solid symbols) and predicted using Eq. for the whole system 
plotted against the number fraction /? = N'^/ {Nj^ + N"^). Different colors 
represent the volume fraction u as shown in the legend. The corresponding 
arrows show K for /3 = 1, i.e. the limit for infinitely small B particles, 
where the measurements are done after removing all the small particles 
and the static assembly consisted of only A. (b) K for the two extreme 
cases: j3 = /3min (solid symbols) and (3 = 1 (empty symbols), both are 
the monodisperse cases with the latter having 5% fewer particles than the 
former. Lines passing through the data is Eq. (25). The dots represents the 
maximum K obtained from (a) for a given density. 
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(a) 


(b) 


Figure 14; (a) Partial bulk modulus (scaled by constant 2{rA)/k) for the 
A A (big symbols) and AB (small symbols) interactions, plotted against the 


number fraction /? = N'^/ (A^J + from the same data as in Fig. 13(a) 
(b) gs calculated using Eq. ( [27l ), where the dashed line is Eq. (I27f with 
constant linear compacity assumption, i.e., x = 1 for the total mixture 
including the rattlers, see Eq. (§. 
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Figure 15: Bulk modulus (scaled by constant 2(r^)/fe) measured using Eq. 
(25) plotted against size ratio tb/ta- The dashed vertical lines represent 
the radio ratio when particle B fills the void by A formed in triangular, 
tetrahedron, square and cubic lattices respectively, as shown in Fig 


A.l 


Different colors represent the volume fraction v as shown in the legend. For 
the two lowest loose states, all B particles are rattlers and hence = in 
the system. 
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Figure A.l: (a) Triangular (b) tetrahedron (c) square and (d) cubic lattices, 
where the small particle of radius r^, rg, and rg respectively is residing 
between bigger particles of radius rA, just touching them. 



C.l.(a) 


C.l.(b) 


Figure C.l: (a) (A)c against ln(j//t'c) using Eq. (C.8). (b) gp measured 
using Eq. (C.IO) with the assumption that the total force on a particle is 
proportional to square of contacts it has with neighbors and is compared 
with the analytical expression of gp in Eq. (C.7). 
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